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We study the metastability and nucleation of the Blume-Capel model on complex networks, in 
which each node can take one of three possible spin variables (—1, 0, 1}. We consider the external 
magnetic held h to be positive, and let the chemical potential A vary between — h and h in a low 
temperature, such that the 1 configuration is stable, and —1 configuration and/or 0 configuration are 
metastable. Combining the heterogeneous mean-held theory with simulations, we show that there 
exist four regions with distinct nucleation scenarios depending on the values of h and A: the system 
undergoes a two-step nucleation process from —1 configuration to 0 conhguration and then to 1 
conhguration (region I); nucleation becomes a one-step process without an intermediate metastable 
conhguration directly from —1 conhguration to 1 conhguration (region 11(1)) or directly from 0 
conhguration to 1 conhguration (region 11(2)) depending on the sign of A; the metastability of the 
system vanishes and nucleation is thus irrelevant (region III). Furthermore, we show that in the 
region I nucleation rates for each step intersect that results in the occurrence of a maximum in the 
total nucleation rate. 

PACS numbers: 89.75.He, 64.60.Q., 05.50.+q 


I. INTRODUCTION 

Complex networks describe not only the pattern discovered ubiquitously in the real world, but also provide a 
unified theoretical framework to understand the inherent complexity in nature 0,0- A central topic in this field is 
to unveil the relationship between the topology of a network and dynamics taking place on it |3|-l£|. In particular, 
phase transitions on complex networks have been a subject of intense research in the field of statistical physics and 
many other disciplines [7J. Extensive research interests have focused on the onset of phase transitions in diverse 
network topologies. Owing to the heterogeneity in degree distribution, phase transitions on complex networks are 
drastically different from those on regular lattices in Euclidean space. For instance, degree heterogeneity can lead to 
a vanishing percolation threshold [8], the whole infection of disease with any small spreading rate [9( , the Ising model 
to be ordered at all temperatures 10-|l2|, the disorder-order transition in voter models [131 ] . synchronization to be 
suppressed mm and different paths towards synchronization in oscillator network 0, spontaneous differentiation 
of nonequilibrium pattern 0, to list just a few. However, there is much less attention paid to the dynamics of a 
phase transition itself on complex networks, such as nucleation in a first-order phase transition. 

Nucleation is a fluctuation-driven process that initiates the decay of a metastable state into a more stable one 
[Tsl . Many important phenomena in nature, like crystallization 0 , glass formation [20], and protein folding 0 , are 
closely related to the nucleation process. In the context of complex networks, the study of the nucleation process is not 
only of theoretical importance for understanding how a first-order phase transition happens in networked systems, but 
also may have potential implications for controlling fluctuation-driven system-wide transitions in real situations, such 
as the transitions between different dynamical attractors in neural networks [22], the genetic switch between high- 
and low-expression states in gene regulatory networks [23l [24j , a new opinion |25| or scientific paradigm formation 
26] as well as language replacement [27|, [28|] in social networks, and spontaneous traffic jamming [29j, synchronization 
HU, cascading failure 0 and recovery [0 close to an explosive phase transition. 

Recently, we have made a tentative step in the study of the nucleation process of the two-state Ising model on 
complex networks, where we have identified nucleation pathways, such as nucleating from nodes with smaller degree 
on heterogeneous networks [34] and a multi-step nucleation process on modular networks 0. In addition, a size-effect 
of the nucleation rate on mean-field-type networks and a nonmonotonic dependences of the nucleation rate on 
the modularity of networks [35] and on the degree heterogeneity |36] were reported. However, many real systems 
possess complicated free-energy landscape with several local minima where phase transition happens usually via these 
intermediate metastable states [0 . The presence of intermediate metastable states has been shown to play a key 
role in determining the pathway and rate of nucleation. For example, it was recently reported that an intermediate 
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metastable phase can provide an easier pathway for the growth of crystal nuclei from fluids, with implications for the 
crystallization of protein and colloid |38l - l42| . Thus, it is natural to generalize the nucleation of the two-state Ising 
model to a three-state spin model on complex networks in which an intermediate metastable state may exist. 

In this paper, we shall use the three-state Blume-Capel (BC) model to investigate the nucleation on complex 
networks. The BC model is a spin-1 Ising model that has been introduced, by Blume and Capel [Hj independently, 
as a model for magnetic systems and then applied to multicomponent fluids [HI]. The BC model defined on a two- 
dimensional lattice has been previously used to study the metastability and nucleation in the limit of zero temperature 
fill and in the absence of external magnetic fieldj47( . Recently, a reentrance phase transition has been observed in the 
BC model defined on heterogeneous networks [48l | . Here, we show, by using mean-field analysis and simulation, that 
there are four distinct regions corresponding to different nucleation scenarios for the networked BC model. Depending 
on the model’s parameters, the system undergoes either a two-step nucleation process with an intermediate metastable 
state or a one-step nucleation process. We also calculate the rates of nucleation by a rare-event sampling method that 
agree with the theoretical predictions by evaluating the free-energy barrier to nucleate. 


II. MODEL 

We consider the BC model defined on a network, where spin variable of each node can take three possible values 
tTi € {—1,0,1}, and interacting according to the Hamiltonian 

n = J Yl M 0 * ~ ~ A a i ~ h J2 ai ’ w 

i<j i i 

where J is the ferromagnetic interaction constant among nodes, A and h have the meaning of the chemical potential 
and the external magnetic field imposed on each node, respectively. The elements of the adjacency matrix of the 
network take ciij = 1 if nodes i and j are connected and a^- = 0 otherwise. 

The present paper is devoted to the study of metastability and nucleation of the networked three-state BC model 
at the low temperature. For the purpose, we first consider the stability of the system in the zero temperature limit. 
In this case, the (local) stable equilibrium refer to the configurations with all the spins equal to — 1, 0, 1, respectively. 
For the sake of convenience, we use —1 , 0, 1 to denote these stable ordered configurations, respectively. Their energy 
are as follows: h — A, 0, and —h — A. Since we want to study the nucleation from — 1 to 0, and then to L we set 
— 1 < 0 < 1 as the relative stabilities of these configurations. To the end, it is required that h — A>0>— h — A, or 
equivalently, —h < A < h and h > 0. Due to the small thermal fluctuation at the low temperature, it is expected that 
the behavior at the low temperature is similar to that at the zero temperature. However, in the presence of small 
thermal fluctuation the notations —1 , 0, 1 refer to the configurations with most instead of all the spins equal to — 1, 
0, 1, respectively. Here, the temperature is fixed at T = 5 (in unit of J/ks) throughout the paper where ks is the 
Boltzmann constant. 


III. THEORY AND SIMULATION 

To proceed the heterogeneous mean-field theory, we first define X as the probability that a node of degree k is in 
the state a € {—1,0,1}. The interaction energy of an edge connecting a fc-degree node and a fc'-degree node is thus 
written as, 


Ekk' — J 


X 


(i) 


( 4 ° ) + 44 " 1) ) 

J ('Qk + Qk’ — 2MfcMfc/) , 
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( 0 ) 
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where 


M fc = x™ - xjr 1] 

Q k = xi 1) +Xi~ 1) 

are the average magnetization and the average squared magnetization of a node of degree k, respectively. 
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Plus the single-node energy, the total energy of the system can be expressed as 

E = P(k)k P(k'\k) E kk ' - A7V ^ P(k)Q k -hN^ P{k)M k 

k k' k k 

= JN (k) (V - M' 2 ) - A NQ - h.NM (4) 

where P(k) is the degree distribution, and P(k'\k) is the conditional probability that a node of degree k links to a 
node of degree k '. We use P(k'\k) = k'P(k')/(k) under the consideration of without degree correlation, where (k) is 
the average degree. In Eq.(4), we have used the definitions, 


]TP(fc)A4 

k 

(5a) 

k 

(5b) 


(5c) 


(5d) 


where M and Q are the average magnetization and the average squared magnetization per node, respectively. M' 
and Q' are the average magnetization and the average squared magnetization of a randomly chosen nearest node, 
respectively. 

Furthermore, let us define S k as the entropy of a node of degree k, the total entropy of the system is 


S = Nj2P(k)S k , 

k 


( 6 ) 


with 


s k = -k B Y J xr in4 • 


(7) 


Combining Eqs.(4) and (6), we can get the expression of free energy, F = E — TS. 
In order to get the extrema of the free energy, we use the Lagrange function, 


= -/+E p ( fc )^ p-Ee 


(a) 

k 


( 8 ) 


where / = F/N is the average free energy per node, and /ifc is the Lagrange multiplier to maintain the normalization 
condition. By minimizing of Eq.(8) with respect to X^\ one has, 
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jj = Jk (1 - 2 M') - A -h 
= 0 


de 

dX^ 

gx t- 1) = '7^ (1 + 2 M') — A + h 


( 10 ) 


where /? = 1 /(fc^T) and e = E/N is the average energy per node. 
Substituting Eq.(9) into Eq.(3), we get 
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2 sinh \/3 (2 JkM' + ft)] 

2 cosh [/3 (2 JkM' + h)\ + exp [/3 (Jk — A)] 
2 cosh \j3 (2 JkM' + h)\ 

2 cosh [/3 (2 JkM' + h)\ + exp [/3 (Jk — A)] 


( 11 ) 

( 12 ) 


Furthermore, inserting Eq.(ll) into Eq.(5c), we get a self-consistent equation of M' that can be numerically solved, 


m> = y; 

k 


kP(k ) 2 sinh [/3 (2 JkM' + h)\ 

( k ) 2 cosh [/3 (2JkM' + h)] + exp [f3 (Jk — A)] 


(13) 


Once the solutions of M' are determined, all quantities will be obtained, including the extrema of free energy. 

To begin with, we consider a Erdds-Renyi (ER) random network with the Poissonian degree distribution P(k ) = 

(fc) fc e - ^ j ft! The average degree we use is ( k ) = 20. By numerically solving Eq.(13), we plot the solutions of 

M' as a function of A for three typical values of ft, as shown in Fig.l. The stable and unstable solutions are depicted 
as the solid and dashed lines, respectively. When h is small, for example ft = 2 in Fig.1(a), there are five solutions 
in the whole allowable range of A € [—ft, ft], in which three of them are stable. Let M' Sl , M' S2 , and M' 3 denote these 
stable solutions around —1, 0, and 1, respectively. These stable solutions give the three stable states —1 , 0, and 1, 
respectively. The other two unstable solutions, M' USi and M ', give the two transition states from ^1 to 0 and from 
0 to 1, respectively. When ft becomes relatively larger, as shown in Fig. 1(b), the stable solution M a and the unstable 
M' USi collide and annihilate each other at A = A C1 via a a saddle-node bifurcation. Meanwhile, M' S2 and the unstable 
M' US2 collide and vanish at A = A C2 > A Cl in the same way. In this case, the property of solutions can be classified into 
three regions depending on the value of A. For A < A Cl , there are two stable solutions, M' S2 and M' 3 , and one unstable 
solution M ' US2 . For A > A C2 , there are also two stable solutions and one unstable solution, but they are M' ai M a , and 
M' us i- While for A Cl < A < A C2 , the property of solutions is the same as in Fig. 1(a). With further increasing ft, as 
shown in Fig. 1(c), A Cl shifts to a larger value, and at the same time A C2 shifts to a smaller value, so that A Cl > A C2 . 
In this case, there is a single stable solution M' at A C2 < A < A Cl . 

For comparison, we have performed Monte Carlo simulations in ER random networks with network size N = 1000 
to compute the steady state values of M'. The simulations start from numerous various initial configurations to 

sample all possible steady state values of M ', where M' can be conveniently computed as M' = ft i&i j ((ft) N) 
with ki being the degree of node i. The simulation results are added into Fig.l as shown by square points. It is clear 
that the simulation results are in excellent agreement with our theoretical estimations. 

To get a global view, we have plotted the phase diagram in the ft ~ A plane, as shown in Fig.2. The phase diagram 
is divided into four different regions according to the property of solutions of M ', separating by the lines of A Cl ~ ft 
and A C2 ~ ft. In region I, M' have five solutions: three of them are stable (corresponding to three stable states 
— 1 , 0, and 1), and the others are unstable (two transition states from ^1 to 0 then to 1). In region 11(1), M' have 
three solutions: two of them are stable, and the other is unstable (two stable states —1 , and 1 and one transition 
state). In region 11(2), M' have also three solutions, but the two stable states are 0, and 1 and one transition state 
between them. In region III, M' has only one stable solution that corresponds to the state 1. Also, we have given 
the simulation results of A Cl (ft) and A C2 (ft), as depicted by square points in Fig.2. One can see that there exist some 
mismatches between the theory and simulation. This is because that near the boundaries the lifetimes of metastable 
states are rather short (or have very low free-energy barrier to nucleate that will be illustrated later), so that such 
metastable states are hard to identify in the simulations. 

We have also performed the calculations in Barabasi-Albert (BA) scale-free networks j2j with the same size and 
average degree, and found that the phase diagram is almost the same as that in ER random network (not shown 
here). That is, the phase diagram is not almost affected by network topology. 

In Fig.3, we have schematically demonstrated the nucleation process in different regions. In region I, the system 
undergoes a two-step nucleation process. The first stage is the nucleating of the nretastale 0 from ^1. Subsequently, 
in the second stage, the transition from 0 to 1 happens via the nucleating of L The free-energy barrier of the two-step 
nucleation are AF_i^o = NAf_i^o = - /m'J and Af^i = N A/o_>i = N(f M ^ S2 - /m'J, respectively. 

Herein, we use the notation /m' to denote free energy per node at M' = M '. In regions 11(1), since the state 0 is no 
longer present, nucleation happens via a one-step process directly from to 1. The resulting free-energy barrier is 
= NAf-i^i = N (/m' ;si — fM ' sl )■ In region 11(2), since the state ceases to exist, nucleation also proceeds 
by one-step process from 0 to I, and the corresponding free-energy barrier is AFo->i- In region III, the only stable 
state is 1 and the nucleation is thus irrelevant. 

Since nucleation rate R is exponentially dependent on — (3AF, R ~ exp(— (3AF), at the top and bottom panels of 
Fig.4 we show that —f3AF as a function of A at three typical different ft: ft = 2, 3.5, 6 (from left to right) in ER 
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FIG. 1: (color online) The solutions of M' as a function of A at three different h: h = 2 (a), h = 3.5 (b), and h = 6 (c). 
The solid and dashed lines depict the stable and unstable solutions of M', respectively, given by the analysis Eq.(15) with a 
Poissonian degree distribution. The square points are given by MC simulations in N = 1000 ER random networks. The other 
parameters are (k) = 20 and T = 5. 


networks (top panels) and BA networks (bottom panels), respectively (shown by lines). To validate the analytical 
results, computer simulation for calculating nucleation rate is desirable. 

However, nucleation is an activated process that occurs extremely slow, and brute-force simulation for observing 
nucleation process is thus prohibitively expensive. To overcome this difficulty, we will employ a recently developed 
simulation method, forward flux sampling (FFS) mm. This method allows us to calculate nucleation rate and 
determine the properties of ensemble toward nucleation pathways. This method uses a series of interfaces in phase 
space between the initial and final states to force the system from the initial state A to the final state B in a ratchet¬ 
like manner. Before the simulation begins, an order parameter r is first defined, such that the system is in state A if 
r < ro and it is in state B if r > r n . A series of nonintersecting interfaces r,; (0 < i < n) lie between states A and 
B , such that any path from A to B must cross each interface without reaching ?y+i before fj. The algorithm first 
runs a long-time simulation which gives an estimate of the flux 4 >a,o escaping from the basin of A and generates a 
collection of configurations corresponding to crossings of interface tq. The next step is to choose a configuration from 
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FIG. 2: (color online) Phase diagram in the h ~ A plane. The phase diagram is divided into four different regions according to 
the property of solutions of M', separating by the lines of A C1 ~ h and A C2 ~ h. In region I, M' have five solutions: three of 
them are stable (corresponding to three stable states —1 , 0, and 1), and the others are unstable (two transition states from 
to 0 then to 1_). In region 11(1), M ' have three solutions: two of them are stable, and the other is unstable (two stable states 
— 1 , and 1_ and one transition state). In region 11(2), M' have also three solutions, but the two stable states are 0, and 1 and 
one transition state between them. In region III, M' have only one stable solution, and the corresponding state is 1. Inset: 
Graphic solution of Eq.(15) for four typical ( h , A) points (stars in phase diagram) chosen in four different regions. Square points 
give the simulation results of X ci {h) and A c 2 (h). The other parameters are the same as those in Fig.l. 


this collection at random and use it to initiate a trial run which is continued until it either reaches rq or returns to 
ro- If rq is reached, store the configuration of the end point of the trial run. Repeat this step, each time choosing a 
random starting configuration from the collection at ro. The fraction of successful trial runs gives an estimate of of 
the probability of reaching r\ without going back into A , P (ri|ro). This process is repeated, step by step, until r n is 
reached, giving the probabilities P (r^+i |fj) (i = 1, • ■ • , n — 1). Finally, we get the nucleation rate R from A to B as 

R = $ A0 P(r n \r 0 ) = ^a.oTT P (r i+ i |a), (14) 

where P (r n |ro) is the probability that a trajectory crossing ro in the direction of B will eventually reach B before 
returning to A. 

For comparison, we have also added the simulation results to Fig.4 (shown by symbols). On one hand, the analytical 
results are in well agreement with the simulation ones. On the other hand, the results in ER random networks and 
in BA scale-free networks are qualitatively the same. 

For h = 2, nucleation is a two-step process in the whole allowable range of A, and the corresponding nucleation rates 
are R- i_>.o and respectively. As A increases, P_i_>.o decreases monotonically, but Rq->i increases monotonically. 

In this case, the total nucleation rate R is expressed as R= {RZi^o + P^i) -1 - It can be seen that R is dominantly 
determined by the smaller of R- i_>o and Therefore, there exists a maximal R at A = A op t where R- i_>.o and 

Ro->i intersect. Here, A op t — 0.1 for ER random networks and A op t — 0.2 for BA scale-free networks that are robust 
to different h. For h = 3.5, nucleation is a two-step process in a certain range of A around zero, while in other range 
nucleation becomes a one-step process. The resulting total nucleation rate has also a maximum at A = A op t- For 
h = 6, nucleation is a one-step process and the corresponding nucleation rate increases as A gradually approaches zero 
until nucleation becomes irrelevant when A crosses the line of X Cl {h) or X C2 (h). 
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FIG. 3: (color online) Schematic plot for nucleation and growth in different regions. In region I, nucleation from —1 to 1 
proceeds in a two-step process via an intermediate metastable state 0. In both regions 11(1) and 11(2), nucleation proceeds in 
a one-step process. In region III, the only state 1 is stable and thus nucleation is irrelevant. 






A A A 


FIG. 4: (color online) The free energy barrier —/3AF and the logarithm of rate of nucleation In 7?, obtained from the HMF 
analysis (lines) and FFS simulations (symbols), as a function of A at three different h\ h = 2, 3.5, 6 (from left to right). Top 
and bottom panels show the results in ER networks and BA networks, respectively. The other parameters are the same as 
those in Fig.l. 


IV. CONCLUSION 

In conclusion, we have studied the nucleation of the three-state BC model on complex networks. By using the 
heterogeneous mean-field theory and simulations, we have found that there exist four distinct regions with different 
nucleation scenarios depending on the model’s parameters: the external field h and the chemical potential A. For 
a small h, or for a moderate h and simultaneously a small |A|, the system goes through a two-step nucleation 
process from configuration to 0 configuration and then to 1 configuration. For a moderate h or a large h and 

















simultaneously a large |A|, nucleation is a one-step process without an intermediate nretastable configuration directly 
from —1 configuration to l configuration or directly from 0 configuration to 1 configuration depending on the sign 
of A. For a large h and simultaneously a small |A|, the nretastability of the system vanishes and nucleation is thus 
irrelevant. Moreover, we have calculated the nucleation rates and found that in the two-step nucleation region there 
exists a maximum for the total nucleation rate. All the results are demonstrated in ER random networks and in BA 
scale-free networks and are qualitatively the same for different network topologies. Quantitatively, the optimal X opt 
at which the total nucleation rate is maximal in ER random networks is less than that in more heterogeneous BA 
scale-free networks. 
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